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ABSTRACT 

We present a new method for optimally extracting point-source time variabil- 
ity information from a series of images. Differential photometry is generally best 
accomplished by subtracting two images separated in time, since this removes 
all constant objects in the field. By removing background sources such as the 
host galaxies of supernovae, such subtractions make possible the measurement 
of the proper flux of point-source objects superimposed on extended sources. In 
traditional difference photometry, a single image is designated as the "template" 
image and subtracted from all other observations. This procedure does not take 
all the available information into account and for sub-optimal template images 
may produce poor results. Given total observations of an object, we show how 
to obtain an estimate of the vector of fluxes from the individual images using the 
antisymmetric matrix of flux differences formed from the N{N — l)/2 distinct 
possible subtractions and provide a prescription for estimating the associated 
uncertainties. We then demonstrate how this method improves results over the 
standard procedure of designating one image as a "template" and differencing 
against only that image. 

Subject headings: methods: data analysis — techniques: photometric — super- 
novae: general 



1. Introduction 



The astronomical time domain provides unique insight into a range of astrophysical 
phenomena. Studies of variable stars yield information about stellar structure and evolution 
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as well as help to set the extra-galactic distance scale. Active Galactic Nuclei (AGN) reveal 
the high-energy phenomena associated with the supcr-massivc black holes that reside at 
the centers of most galaxies. Supernovae (SNe) and Gamma-Ray Bursts (GRBs) provide a 
glimpse of the fantastic energies released during the violent death throes of several types of 
stars. Type la supernovae (SNe la) are of particular interest because their use as "standard 
candles" has revealed the acceleration of the expansion of the universe from an inferred 
cosmological constant-like force (Riess et al. 1998; Perlmutter et al. 1999). 

Studies of variable sources require specific analysis methods that are not necessary for 
non-variable sources. Since it is often difficult to detect variation in an object by simply 
inspecting images, the standard procedure is to subtract images taken at different times to 
remove objects with constant flux. Photometrically variable objects are then obvious. For 
the case of SNe, one typically obtains a pair of observations separated in time to allow for 
SNe not present in the first image to reach observable brightness in the second (see, e.g., 
Perlmutter et al. 1995 and Schmidt et al. 1998 for a description of the method). After 
detection, additional observations are made to obtain the complete light curve necessary 
for cosmological analysis (see Phillips 1993; Riess, Press, & Kirshner 1996). In order to 
construct the light curve, it is necessary for at least one observation (the "template" image) 
to contain no SN flux. Often this is the initial image used during discovery of the SN. In 
many instances, however, the SN is present at a faint level in this image, so an additional 
observation, taken after the SN has faded from view, is required. The fight curve is then 
calculated by measuring the flux levels in subtractions of each image from the designated 
template using, for example, the subtraction procedure described by Alard & Lupton (1998). 
This, the "single-template method," is the typical means of constructing light curves of SNe 
and other variable sources. 

However, this method has certain drawbacks. The primary flaw is that the quality of 
any subtraction depends greatly upon the two images involved. If the template is of a poor 
quality caused, for instance, by poor seeing or a low signal-to-noise ratio (S/N), then every 
subtraction will be degraded, with a corresponding increase in the measured flux uncertainty, 
even if all other images are of high quality. Any flaw in the template creates a systematic 
error for the entire light curve that is not detectable from internal consistency checks or 
through comparison with another SN light curve. 

In order to alleviate this problem, we have developed a new method for constructing 
light curves of photometrically variable objects. Given N observations there arc a total of 
A^(A^— 1)/2 pairs of images that can be subtracted together, only A^— 1 of which are performed 
in the single-template method. A matrix of flux differences can be constructed from these 
subtractions and used to determine the flux at each individual epoch. This process removes 
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the dependence on any single observation, because all observations are treated equally as 
a "template." We refer to this method as the "N(N-l)/2" method (hereafter abbreviated 
NN2; see Novicki & Tonry 2000 for an initial description). 

Section 2 describes the mathematical underpinnings of NN2. In Section 3 we demon- 
strate the efficacy of the method using simulated SNe inserted into images used during an 
actual high-redshift SN survey. Section 4 gives our conclusions. 

2. Mathematical Basis of the NN2 Method 

We assume that we start with N observations of an object, so that one may construct 
from all pairs of subtractions an N x N antisymmetric matrix A of flux differences that we 
wish to analyze as a "vector-term difference." In other words, we want to flnd an A'"-vector 
V of fluxes such that 

= Vj - V,. (1) 

We also assume that we have a symmetric N x N error matrix E that expresses our uncer- 
tainty in each term of A. As we shall see, this matrix may not be easy to generate, and its 
interpretation may be somewhat ambiguous. However, one can imagine generating an error 
matrix by the following procedure. 

In each of the difference images, we measure a flux for the object in question. In 
general this measurement consists of fitting a fixed point-spread function (PSF) profile at 
the location of the object by adjusting the amplitude of the PSF (both positive or negative) 
and the local background level. The PSF profile may be obtained from a suitable star in 
the original image while the location of the variable source may be determined by summing 
all the difference images (adjusted to keep the sign of the object positive) and fitting the 
location in this sum. Once we have a flux measurement, we can insert copies of the object 
at nearby empty regions of the difference image and repeat the procedure. The mean of the 
recovered fluxes indicates whether there is a bias in the measurement, and the scatter may 
be used as a term Eij in the error matrix. 

The crux of the NN2 method is the distillation of the photometric measurements from 
the full set of N{N — l)/2 subtractions to a lightcurve, V, that represents our best under- 
standing of the behavior of the object under consideration. As long as it is consistent, the 
exact procedure for measuring the flux on the difference images is not central to the NN2 
method we present here. 
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In order to find an optimal V, we wish to minimize the quantity 



(2) 



This construction may not be entirely appropriate depending on the errors in the flux dif- 
ferences Aij. Ideally, if we possessed an extremely high-quality template with no SN flux 
present and applied an optimal subtraction procedure, the errors would be primarily due 
to photon counting statistics (see Alard & Lupton 1998 for a discussion). We would ex- 
pect these errors to be uncorrelated and would simply wish to employ the single-template 
method to construct the SN light curve. As mentioned in the introduction, in practice there 
are nearly always imperfections associated with the template image that remove us from 
this idealized regime. These template errors introduce correlations in the individual flux 
measurement errors that are difficult to quantify and are typically assumed to be negligible 
in SN light-curve analysis. The use of the NN2 method, however, will introduce further 
correlations as a result of the common images in the various subtractions (for instance, the 
error in Vx — V2 will be anti-correlated with the error in V2 — V3 due to the common error 
in V2). Although we believe that the use of the NN2 method will improve the ability to 
accurately recover variable object light curves, one should recognize that the NN2 method is 
expected to introduce these additional correlations to the fluxes measured from the various 
subtraction images, and so the given above is not technically appropriate. Errors due to 
systematics in the subtraction procedure, such as those associated with template or software 
imperfections, would be expected to be effectively uncorrelated, and if they were dominant 
then Eq. 2 would indeed represent the proper x^- With these caveats in mind, we will pro- 
ceed to use the definition of x"^ as given in Eq. 2 as the basis of the NN2 method. Tests of 
its ability to recover accurate light-curve information in the following section demonstrate 
its effectiveness in practice. 

However, we need to make one minor modification to our x^ because the defined 
in Eq. 2 is degenerate to the addition of a constant to the V vector — geometrically, is 
constant along the line In order to lift this degeneracy and permit us to solve for V, 

we add a term to x^ thai is quadratic in the degenerate direction, so that 




(3) 



where (E) is a suitable typical uncertainty; for example. 




(4) 
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Our solution will therefore have Vi — 0. This construction exphcitly forces one to deter- 
mine an accurate zero flux level at a later stage. In the single-template method this zero 
flux level is generally implicitly determined by assuming that the object of interest has zero 
flux in the template image. This same assumption can similarly be used in the NN2 method, 
but more sophisticated methods involving comparisons of many different images can also 
be invoked. If the absolute brightness of the variation being studied is important, the NN2 
method clearly cannot free one from the requirement of having a fiducial image to measure 
the zero flux level. This is a fundamental limitation of any differential photometry method as 
the information is simply not available without such a fiducial image. However, even in the 
absence of a fiducial image, the NN2 method will produce a sensible and accurate relative 
lightcurve. 

We now seek to solve for our lightcurve vector V by minimizing with respect to V: 

i,j;i<j i ^ ' 

Exploiting the antisymmetry of A and the symmetry of E we can rewrite Eq. 6 as 
These equations can be solved for V by inverting a matrix C: 

where 111 

The inverse of this curvature (Hessian) matrix C yields uncertainties for V from the 
square root of the diagonal elements as well as covariances from normalizing the off-diagonal 
elements by the two diagonal terms (under the assumption that the error matrix truly does 
represent Gaussian, independent uncertainties for each of the terms of the antisymmetric 
difference matrix). 

An alternative approach to calculating uncertainties in V stems from assuming that 
there is a vector a such that 

Ei^af + a]. (10) 
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Under this assumption, we seek to minimize 

XI=E {-Ei + a^ + a])\ (11) 

i,j;i<j 



The minimization condition is 



= ll (12) 



'k 

= 2 J2 {-Ei + 'yI + 'y]){S^k + SJk) (13) 
- 2 5^(-£;^ + a,^ + a^). (14) 

These N equations are solved by inverting a matrix D 
where 

Afe = 1 + (iV - 2)5,fc. (16) 

After solving for V and cr, we can evaluate the quality of the fit by comparing to the 
number of degrees of freedom, 

This N^of comes from the number of data points, N{N — l)/2, minus the number of model 
parameters, N — 1. Recall that we've explicitly required '^^Vi = 0, so that the number of 
model parameters is A^" — 1 rather than N. 

Having outlined the basic method, wc now discuss a fundamental uncertainty in the NN2 
process. We can imagine two types of error that will cause V to differ from the true flux 
values. The first, which we term "external error," is intrinsic to the images themselves. For 
example, if the object has a positive statistical fiuctuation in fiux in one image or is corrupted 
by a cosmic ray that happens to coincide with the position of the object on the detector, this 
error will propagate through the entire differencing and analysis procedure. It is possible to 
obtain an antisymmetric difference matrix that is an exact vector-term difference (x^ = 0), 
but the solution vector will still contain errors. The second type of error, which we call 
"internal error," is caused by the procedure of generating the antisymmetric matrix. One 
might imagine a set of images that have no ffux error whatsoever, but through errors in 
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convolving, differencing, or fiux fitting, an antisymmetric matrix may be created tfiat is not 
a perfect vector-term difference and for whicfi > 0. 

Rougfily speaking, one might expect that if the error matrix E consists entirely of ex- 
ternal errors the resulting a terms will all be approximately E/y/2, since E is the quadrature 
sum of two a terms. Alternatively, if the error matrix is purely internal error the a terms 
might be expected to be approximately E /\fN , since each term in V comes from comparison 
with A^ — 1 other images. In the case of external errors, the uncertainties derived from the 
analysis are correct. In the internal error case the uncertainties obtained from the covariance 
matrix derived from the analysis are likewise appropriate. 

It is not clear how to disentangle these different sorts of errors. The procedure suggested 
above of dropping copies of the object into each difference image and evaluating the scatter 
of the result will be sensitive to each sort of error, but it is possible to imagine cases where 
this procedure is unsatisfactory. We suggest that the errors provided in the E matrix be 
interpreted as external errors and taken seriously as such. Thus, the vector V is assigned 
an external uncertainty equal to the cr vector. However, in order to handle a situation 
where x^/^dof is much greater than 1 (i.e., where the antisymmetric matrix is simply not 
well represented as a vector-term difference) , we suggest also creating an internal uncertainty 
vector T that is obtained from the diagonal terms of the covariance matrix, scaled by X^/^dof' 

The total uncertainty is then the quadrature sum of a and r. Note that this approach 
implicitly assumes that the internal and external errors are uncorrelated and are proportional 
to one another as well as the provided E matrix. 

For problems where X^/^dof is near unity without adjustment, the r vector will be 
smaller than the a vector by approximately \/2/N and will make a fairly small contribution 
to the total uncertainty. When x'^/^dof ^ 1 (i-e., the antisymmetric matrix is very closely 
represented by the vector-term difference), the r vector will be negligible. However, when 
X^/Nfiof ^ Ij the T vector will act to correct x^/iVdo/ to approximately unity, and this 
procedure will provide reasonable uncertainties, even though E may be much too small. 
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3. Demonstration of Improved Accuracy in Recovering Variable Object Light 

Curves 

The first extensive use of tlie NN2 metliod we fiave developed liere occurred during 
tlie SN-searcli component of tlie If A Deep Survey (Barris et al. 2004), altliougli we also 
employed it to a limited extent in a previous SN survey by Tonry et al. (2003). The 
IfA Deep Survey was undertaken primarily with Suprime-Cam (Miyazaki et al. 1998) on 
the Subaru 8.2-m telescope and was supplemented with the 12K camera (Cuillandre et al. 
1999) on the Canada-Prance-Hawaii 3.6-m telescope. Scores of high-redshift SN candidates 
were discovered (Barris et al. 2001, 2002) with 23 confirmed as SNe la. We here present 
several tests we performed to demonstrate the improved performance of NN2 relative to the 
single-template method. 

In order to make a controlled test of the effectiveness of NN2 vs. the single-template 
method, we inserted artificial SNe into the survey images. The light curves of these objects 
consisted of a linear ramp-up and ramp-down in brightness over the time period covered by 
the survey observations. The timing of the light-curve maxima were selected at random and 
could lie within or outside of the survey period. The simulated SNe were laid down in a 
regular grid across the survey area, and all pairs of images were subtracted. Object-detection 
software was then run on all subtraction images to detect photometrically variable objects 
(both real objects and the artificial SNe) and to construct the NN2 flux difference matrix. 
The positions of both real and artificial SNe were fit by the object-detection software as 
described in Section 2. Since we knew the true light-curve properties used to create the 
synthetic SNe, we could calculate the root-mean-square scatter (RMS) around this artificial 
light curve using both the NN2 fiux calculation and the single-template method with every 
individual observation as the template (this latter is equivalent to taking the fiux values from 
a single column of the NN2 flux difference matrix) . 

We inserted approximately 2000 simulated SNe into the J-band observations of each of 
four ~0.5 square-degree fields from the IfA Deep Survey, spanning a peak magnitude range 
of approximately rrij = 21-25. We used a predefined grid of positions to insert the simulated 
SNe, without taking into consideration the presence of actual objects nearby that would 
cause problems for detection and accurate photometric measurements. The small fraction 
that were so affected were accordingly not used in the final analysis. 

Figure 1 shows the percentage improvement in the cumulative distributions of the RMS 
(in flux units) from the NN2 method over the set of all RMS values from the single-template 
method from the four survey flelds (RMS values are calculated from flux measurements scaled 
so that a value of flux = 1 corresponds to a magnitude of 25). The cumulative fraction for 
the NN2 method is larger than that for the single-template method distribution at all values 
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of RMS, indicating that the NN2 method does indeed tend to yield smaller RMS values. 
The NN2 method more accurately recovers the actual hght curve of these variable objects. 

Since one could imagine that certain templates of very high quality could outperform 
the NN2 method while the collection of all single-template measurements, as shown in Fig. 1, 
does not, we next examine the relationship between the NN2 RMS and the single-template 
RMS for individual observations. To illustrate this comparison, we will concentrate on only 

one of the survey fields, although details of our investigation of the entire survey area can be 
found in Barris (2004). Table 1 contains relevant information for the 16 observations of the 
selected field, f0438. This field is representative of the entire survey area, though it is notable 
that it contains an observation that was quite strongly affected by clouds (Observation 4), 
as seen by its unusually bright zero-point magnitude. Also noteworthy is Observation 11, 
taken in poor seeing conditions. We would expect the performance of the single-template 
method using these observations to be poor in comparison to the NN2 procedure. 

In Table 1 we demonstrate that the typical RMS obtained with the NN2 method for 
the set of 1775 simulated SNe is smaller than the single-template method RMS using every 
observation of the selected field. The improvement is generally fairly small, ranging from 
~ 5 — 10%. We demonstrate below that these differences are statistically significant. The 
use of either Observations 4 and 11 as single templates, as expected, produces substantially 
worse results relative to the NN2 method than the other observations. For these observa- 
tions the improvement due to the NN2 method is substantially larger than 10%. Figure 2 
shows graphically the percentage difference in the cumulative RMS distributions, similar to 
Figure 1, for each individual observation of f0438. 

Having demonstrated the improved performance of the NN2 method, we can test the 
statistical significance of the differences between the NN2 RMS values and those calculated 
via the single-template method and examine whether these differences indicate an actual 
difference in the distributions of the results from the two methods. To do so we use the non- 
parametric Kolmogorov-Smirnov (K-S) test, with results given in Table 2. For the sample of 
all single-template RMS values compared to NN2, the K-S probabihty value is < 5 x 10~^, 
indicating with strong confidence that the distributions are different. We also divide the 
sample into magnitude bins, since the relative behavior of NN2 RMS to single-template 
RMS is expected to be sensitive to the object's S/N and hence to the magnitude for a given 
sensitivity. 

The K-S probability values for three approximately equal magnitude bins show that 
for each of the subsamples the difference between NN2 and the single-template method 
is statistically significant, increasingly so at fainter magnitudes. Finally, we compare each 
observation individually with the NN2 method and see again that the observed improvement 
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in RMS with NN2 is highly significant for nearly all observations (the only obvious potential 
exception is observation 8). These K-S probability values demonstrate that the NN2 method 
is not distributed identically to the single-template method, and the differences in median 
values given in Table 1 are indeed indicative of statistically significant differences in the 
distributions. This test confirms that the NN2 method truly does produce improved results 
in generating differential light curves. 

4. Conclusions 

We have described the mathematical foundation of a new method for constructing the 
light curves of photometrically variable objects. This method uses all N{N — l)/2 possible 
subtractions involving images in order to calculate a vector of fiuxcs of the variable object 
and offers a powerful alternative to the single-template method that is in standard use for 
studying variable sources. 

If one has a data set with a limited number of good fiducial observations, the NN2 
method will outperform any single-template subtraction approach. For cases where a large 
number of fiducial observations are available to construct a deep template image, the NN2 
method and the single-template approach using this deep template should yield comparable 
results. In this situation we would encourage the use of both methods to provide additional 
checks and constraints on the differential light curve. 

We have tested the performance of the NN2 method by inserting artificial SNe into 
images from the IfA Deep Survey and comparing the RMS scatter from flux measurements 
using the two different methods. We flnd that the RMS from the NN2 method is better 
than the single-template RMS for the large majority (typically 65%-72%) of the SNe for 

every possible template. The median values for the ratio of NN2 RMS (in fiux units) to 
single-template RMS measurements arc typically 0.93-0.96, demonstrating that the NN2 
method results in a ~ 5% improvement in the accuracy of the recovered light curve for these 
observations. 

Using Kolmogorov-Smirnov statistics, we have demonstrated that these differences are 
significant, refiecting an actual difference between the performance of the two methods. We 
find extremely high probabilities that the NN2 RMS is distributed significantly differently 
from the single-template RMS values. This difference and improvement in RMS holds even 
for the very high quality templates that would be considered ideal for the single-template 
method. 

We therefore make the following conclusions: 
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1. For the If A Deep Survey observations, use of the NN2 method typically results in 
a 5-10% improvement in the RMS of the recovered light curve in comparison to the single- 
template method. 

2. For observations that have a large external error, such as those taken under poor 
conditions, the NN2 method results in a substantial improvement 10%) over the single- 
template method. 

3. When working with high-quality observations, with small external error, the internal 
errors (such as those due to implementation of the subtraction process) dominate. If these 
errors are large, the NN2 method should outperform the single-template method to a large 
degree. If these errors are kept small, as we believe is possible based our extensive experience 
with SN surveys, then the NN2 method will result in a modest but significant improvement 
in accuracy of light-curve recovery. 

In summary, the NN2 method we present here maximizes the time variability informa- 
tion contained in a series of observations by using the relative differences between all pairs 
of images to construct the optimal differential light curve, references 

The source code for our implementation of the NN2 method presented here is available 
at http:/ /www.ctio. noao.edu/essence/nn2/. 

This work was supported in part by grant AST-0443378 from the United States National 
Science Foundation. 
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Observed RMS Flux Error (in=25) 



Fig. 1. — The improvement in the cumulative RMS distribution for the NN2 technique 
compared with that of the single-template method for each of four survey fields is shown 
here as the extra fraction of NN2-derived light curves, compared to the conventional light 
curves, that have an observed RMS flux error less than the abscissa value. For each field the 
NN2 cumulative fraction is larger at all RMS values, demonstrating its improved performance 
in accurate light-curve recovery compared to the single-template method. 
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Fig. 2. — The percentage improvement in the cumulative RMS distribution for the NN2 
technique compared with that of the single-template method for each of the 16 observations 
of a representative field, f0438, is shown here as in Fig. 1. Particularly noteworthy are 
the poor results of using as the single template either observations 4 or 11, which suffered 
from heavy cloud cover and poor seeing, respectively. The vertical scale of the sub-plot for 
Observation 4 has been accordingly scaled up by a factor of 3. For templates of higher quality, 
the difference in relative performance is less significant. Nonetheless, NN2 outperforms the 
single-template method for any choice of the single template image. 
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Table 1. F0438 Observational Information and Light-Curve Recovery Comparison 



Tin 


MJD 




zp2 


Median RMS^ 


Median('RMS^r^n/RMS^ 1^ 


1 


52164.57 


0.75 


30.48 


0.441 


0.951 


2 


52176.49 


0.72 


30.05 


0.432 


0.963 


3 


52191.50 


0.71 


30.12 


0.427 


0.970 


4 


52198.56 


0.59 


27.63 


0.810 


0.625 


5 


52204.59 


0.63 


30.46 


0.424 


0.973 


6 


52225.39 


0.62 


30.17 


0.445 


0.960 


7 


52231.38 


0.51 


30.46 


0.455 


0.944 


8 


52232.35 


0.68 


30.21 


0.419 


0.970 


9 


52236.34 


0.94 


30.36 


0.440 


0.938 


10 


52252.35 


0.83 


30.59 


0.429 


0.957 


11 


52263.38 


1.14 


30.33 


0.515 


0.839 


12 


52283.39 


0.68 


30.51 


0.460 


0.946 


13 


52288.23 


0.84 


30.65 


0.438 


0.941 


14 


52289.39 


0.93 


30.27 


0.448 


0.936 


15 


52323.35 


0.97 


30.23 


0.468 


0.908 


16 


52369.25 


0.58 


29.58 


0.432 


0.943 



^Arcseconds 

^Magnitude zero-point of the processed images from each observation 

^RMS calculated for flux units, with flux=l corresponding to m = 25. Median RMS 
value for the NN2 method is 0.422. 

^Median value of the ratio of the RMS as calculated with the single-template method 
to the RMS as calculated with the NN2 method. A total of 1775 simulated SNe were 
inserted in the survey images for this field. 
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Table 2. Field 0438 Kolmogorov-Smirnov Statistics 



Nl 


N2 


cut criterion 


K-b statistic 


K-b probability^ 


1 77^ 


984(1(1 


all 

cm 


076Q 


4 7430e-09 


482 


771 2 


rr? < 23 


0809 


4 9083P-03 


fil 7 


Q872 


23 < rr? < 24 5 


08Q7 


1 ^707p-04 


676 




24 5 < rr? 


1330 


2 5150P-10 


1 775 


1 775 




0524 


1 4698p-02 


1775 


1775 


nli*?p"rvpif inn 9 


0462 


4 3869e-02 


1775 


1775 


observation 3 


0.0411 


9.6888e-02 


1775 


1775 


observation 4 


0.3977 


<1.0e-30 


1775 


1775 


observation 5 


0.0445 


5.7717e-02 


1775 


1775 


observation 6 


0.0608 


2.6510e-03 


1775 


1775 


observation 7 


0.0845 


5.6244e-06 


1775 


1775 


observation 8 


0.0332 


2.7611e-01 


1775 


1775 


observation 9 


0.0518 


1.6326e-02 


1775 


1775 


observation 10 


0.0417 


8.9127e-02 


1775 


1775 


observation 11 


0.1639 


2.5672e-21 


1775 


1775 


observation 12 


0.0794 


2.4901e-05 


1775 


1775 


observation 13 


0.0575 


5.4234e-03 


1775 


1775 


observation 14 


0.0710 


2.4226e-04 


1775 


1775 


observation 15 


0.0885 


1.6592e-06 


1775 


1775 


observation 16 


0.0541 


1.0650e-02 



^The probability that identically distributed distributions will ex- 
hibit a value for the K-S statistic larger than that observed. 



